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THE  FOLLOWING  THREE  PATTERNS  OCCUR  FREQUENTLY. 

BR=WR*AR-WI*AI 

BI=AI*WR+AR*WI 

DATA ( J) =DATA ( I ) -TEMPR 
DATA ( J+l ) =DATA ( 1+1 ) -TEMPI 
DATA ( I ) =DATA ( I ) +TEMPR 
DATA ( 1+1 ) =DATA ( 1  +  1 ) +TEMPI 

INDEX2MAX=INDEX1+N1-N2 

P.  15*  L.  7 
7  ISTEP=2*MMAX 

P.  21 ►  L.  2  AND  P.  17*  L.  2 
2  NTOT=NTOT*NN(IDIM) 

P.  22 *  L.  5-2  AND  P.  17*  L.  100-2 
NP2=NP1*N 

P.  22 *  L.  12  AND  L.  51 
12  OR  51  NTWO=NTWO+NTWO 

P.  22 *  L.  70+2 
I1RNG=NP1 

IF ( IDIM-4 >71*100*100 

P.  23 *  L.  72+1 
I 1RNG=NP0* ( l+NPREV/2 ) 

P.  23*  L.  120  AND  P.  17*  L.110 
110  OR  120  IlMAX=I2+NPl-2 

P.  23*  L.  120+3  AND  P.  17*  L.  110+3 
J3=J+I3-I2 

P.  23*  L.  200 
200  NW0RK=2*N 

P.  23*  L. 210-1 

IF ( I CASE-3) 210  * 220  *  210 

P.  23*  L.  240+1 
J=J+IFP1 

IF ( J-I3-IFP2 ) 260  *  250  *250 

P.  24*  L.  420+1  AND  P.  18*  L.  420+1 
KMIN=IPAR*M+I1 

P.  24*  L.  440  AND  P.  18*  L.  440 
440  KDIF=IPAR*MMAX 
450  KSTEP=4*KDIF 

P.  24*  L.  520+1  AND  P.  18*  L.  520+1 

KMIN=4*  (KMIN-ID  +  Il 

KDIF=KSTEP 

IF ( KDIF-NP2HF) 450  * 450  *  530 

P.  25*  L.  550+1  AND  P.  19*  L.  550+1 
WR=(WR+WI)*RTHLF 
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P.  25f  L.  560+2  AND  P.  19f  L.  560+2 
WI  = (TEMPR+WI ) *RTHLF 

P.  25 '  L.  570+2  AND  P.  19r  L.  570+2 

mmax=mmax+mmax 

P.  26 '  L.  650+2 
J2RNG=IFPl*(l+IFACT(IF)/2) 

P.  26 f  L.  655-2 
I=1+(J3-I3)/NP1HF 

P.  26 f  L.  665 

665  ICONJ=l+ ( IFP2-2+J2+I3+J3) /NP1HF 

P.  27 f  L.  670  +  1 
TEMPI=SUMI 

SUMR=TWOWR*SUMR-OLDSR+DATA(J) 

SUMI=TWOWR*SUMI-OLDSI+DATA(J+l) 

OLDSRrTEMPR 

OLDSI=TEMPI 

J=J-IFP1 

IF ( J-JMIN) 675 f  675  f  670 
675  TEMPR=WR*SUMR-OLDSR+DATA(J) 
TEMPI=WI+SUMI 
WORK ( I ) =TEMPR-TEMPI 
WORK ( ICONJ) =TEMPR+TEMPI 
TEMPR=WR*SUMI-OLDSI+DATA ( J+l ) 
TEMPI=WI*SUMR 
WORK ( 1+1 ) =TEMPR+TEMPI 
WORK ( ICONJ+1 ) =TEMPR-TEMPI 

P.  27 1  L.  690+2 
I2MAX=I3+NP2-NP1 

P.  27 f  L.  710-2 
JMIN=2*NHALF-1 

P.  28 f  L.  740 
740  NP2=NP2+NP2 

P.  28f  L.  745-1 
IMAX=NT0T/2+l 
745  IMIN=IMAX-2*NHALF 

P.  28 f  L.  805+1 
I2MAX=I3+NP2-NP1 

P.  28f  L.  805+3 
IMIN=I2+I1RNG 
IMAX=I2+NPl-2 
JMAX=2*I3+NP1-IMIN 

P.  28 f  L.  810 
810  JMAX=JMAX+NP2 
820  IF ( I DIM-2 )850f850f830 

830  Jr JMAX+NPO 

P.  28 f  L.  840 
840  JrJ-2 

P.  28f  L.  860 
860  JrJ-NPO 
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This  note  describes  and  lists  three  programs,  all  written  in  USASI  Basic 
Fortran,  which  perform  the  discrete  Fourier  transform  upon  a  multi-dimensional 
array  of  floating  point  data.  The  data  may  be  either  real  or  complex,  with  a 
savings  in  running  time  for  real  over  complex  (see  Timing).  The  transform 
values  are  always  complex  and  are  returned  in  the  array  used  to  carry  the 
original  data.  The  running  time  is  much  shorter  than  that  of  any  program 
performing  a  direct  summation,  even  when  sine  and  cosine  values  are  precalcu¬ 
lated  and  stored  in  a  table.  For  example,  on  a  CDC  3300  with  floating  point 
add  time  of  six  microseconds,  a  complex  array  of  size  80  x  80  can  be  trans¬ 
formed  in  19.2  seconds.  Besides  the  main  array,  only  a  working  storage  array 
of  size  160  need  be  supplied. 

The  exact  operation  performed  is  called  finite  discrete  Fourier  trans¬ 
formation,  also  known  as  harmonic  analysis  or  trigonometric  interpolation. 
Given  an  array  of  data  DATA(ll,I2, . . . ) , 

TRMSF0BM(J1,J2,...)  =  2:  [DATA(I1,I2, . . . ) 

W2(I2-1)(J2-1)_}  ( 

where  W1  =  exp(-2iti/Nl) ,  W2  =  exp(-2rti/N2) , . . .  and  II  and  J1  run  from  1  to 
Nl,  12  and  J2  run  from  1  to  N2,  etc.  The  Fortran  convention  of  subscripts 
beginning  at  one  is  adhered  to.  This  summation  possesses  many  of  the  proper¬ 
ties  of  the  more  usual  infinite  integral 

F(y)  =  /  f(x)  e"2:rtlxy  dx  . 

—OO 

By  interpreting  the  subscripts  modulo  Nl,  N2,  etc.  and  requiring  the  data  to 
represent  equispaced  points,  we  can  easily  prove  the  usual  properties  about 
linearity,  orthogonality,  inverse  transform  and  relationship  to  convolution. 
See  Gentleman  and  Sande  ([3]?  1966). 
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There  is  no  limit  on  the  dimensionality  (number  of  subscripts)  of  the  data 
array.  A  three-dimensional  transform  can  be  performed  as  easily  as  a  one¬ 
dimensional  transform,  though  in  a  proportionately  greater  time.  An  inverse 
transform  can  be  performed,  in  which  the  sign  in  the  exponentials  is  +, 
instead  of  -  .  If  an  inverse  transform  is  performed  upon  an  array  of  trans¬ 
formed  data,  the  original  data  will  reappear  multiplied  by  W1*N2*. . .  . 

The  length  of  each  dimension  may  be  any  integer,  and  as  large  as  storage 
will  permit.  However,  the  program  runs  faster  on  composite  integers  than  on 
primes,  and  is  particularly  fast  on  numbers  rich  in  factors  of  two.  For 
example,  on  the  CDC  3300,  the  following  timings  for  a  one -dimensional  transform 
have  been  calculated  from  the  timing  formula: 


I 

Factorization 

Time  for  Complex  Transform  (sec) 

4094 

2  X  23  X  89 

80 

4095 

32  x  5  x  7  x  13 

24 

4096 

212 

6.2 

4097 

17  X  24i 

180 

4098 

2  x  3  X  683 

480 

4099 

prime 

2868 

4100 

22  X  52  X  41 

39 

Calling  Sequence 

The  listings  of  three  programs  are  given  in  the  appendices.  F0UR1  is  a 
subset  of  F0UR2,  which  in  turn  is  a  subset  of  FOURT.  FOURT  is  the  most  general, 
accepting  multidimensional  arrays  of  any  size.  F0UR2  is  the  same  speed  as 
FOURT  but  accepts  only  complex  multidimensional  arrays  whose  dimensions  are 
powers  of  two.  F0UR1  is  much  slower  than  FOURT  or  F0UR2,  and  performs  only 
one-dimensional  transforms  on  complex  arrays  whose  lengths  are  powers  of  two. 
F0UR1  is  intended  mainly  for  pedagogical  purposes;  it  is  half  a  page  of 
Fortran,  the  others  being  much  longer. 
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The  calling  sequences  are: 

CALL  FOURT  ( DATA ,  NN ,  NDIM ,  IS IGN , IFORM , WORK) 

CALL  F0UR2  ( DATA, NN, NDIM, ISIGN) 

CALL  FOUR1  (DATA,NN, ISIGN) 

In  all  cases,  DATA  is  the  array  used  to  hold  the  real  and  imaginary  parts 
of  the  input  data  and  the  transform  values  on  output.  The  real  and  imaginary 
parts  of  a  datum  must  be  placed  into  immediately  adjacent  locations  in  storage. 
This  is  the  form  of  storage  used  by  Fortran  IV,  and  may  be  accomplished  in 
Fortran  II  by  making  the  first  dimension  of  DATA  of  length  two,  referring  to 
the  real  and  imaginary  parts.  If  the  data  placed  in  DATA  on  input  are  real, 
they  must  have  imaginary  parts  of  zero  appended.  The  transform  values  are 
always  complex  and  replace  the  input  data.  Hence,  the  array  DATA  must  always 
be  of  complex  format . 

For  F0TJR1,  array  DATA  must  be  one-dimensional,  of  length  NN.  For  F0UR2 
and  FOURT,  it  may  be  multidimensional.  The  extent  of  each  dimension  (except 
for  the  possible  first  dimension  referring  to  the  real  and  imaginary  parts) 
is  given  in  the  integer  array  NN,  which  is  of  length  NDIM,  the  number  of 
dimensions.  That  is,  NN(l)  =  Nl,  NN(2)  =  N2,  etc.  * 

ISIGN  is  an  integer  used  to  indicate  the  direction  of  the  transform.  It 
is  minus  one  to  indicate  a  forward  transform  (exponential  sign  is  -)  and  plus 
one  to  indicate  an  inverse  transform  (sign  is  +).  The  scale  factor  l/(Nl*N2*. . . ) 
frequently  seen  in  definitions  of  the  Fourier  transform  must  be  applied  by 
the  user. 

If  the  data  being  passed  to  FOURT  are  real  (i.e.,  have  zero  imaginary 
parts),  the  integer  IFORM  should  be  set  to  zero.  This  will  speed  execution 
(see  Timing).  For  complex  data,  IFORM  must  be  plus  one. 

WORK  is  an  array  used  by  FOURT  when  any  of  the  dimensions  of  DATA  is  not 
a  power  of  two.  Since  F0UR2  and  F0UR1  are  restricted  to  powers  of  two,  WORK  is 
not  needed.  If  the  dimensions  of  DATA  are  all  powers  of  two  in  FOURT,  WORK 
may  be  replaced  by  a  zero  in  the  calling  sequence.  Otherwise,  it  must  be 


As  usual,  the  first  subscript  varies  the  fastest  in  storage  order. 
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supplied,  a  real  floating  point  array  of  length  twice  the  longest  dimension 
of  DATA  which  is  not  a  power  of  two.  In  one  dimension,  for  the  length  not  a 
power  of  two,  WORK  occupies  as  many  storage  locations  as  DATA.  If  given,  it 
may  not  be  the  same  array  as  DATA. 

Double  precision  versions  of  these  programs  may  be  obtained  by  changing 
the  names  to  DFOURT,  DF0UR2,  and  DF0UR1,  declaring  double  precision  all 
variables  not  beginning  with  the  letters  I,  J,  K,  L,  M  or  N,  changing  the 
references  to  COS  and  SIN  to  DCOS  and  DSIN  and  assigning  the  correct  precision 
constants  to  TWOPI  (2rt)  and  RTHLF  (0.52)*  DATA  and  WORK  must  then  be  double 
precision  arrays. 

Storage  and  Common 

No  common  of  any  kind  is  used.  An  integer  array  of  length  thirty- two  is 
used  by  FOURT.  FOURT  is  about  four  hundred  Fortran  statements  long,  F0UR2 
about  one  hundred  and  twenty  and  F0UR1  thirty- seven. 

Return  and  Error  Messages 

There  are  no  error  messages,  error  halts  or  error  returns  in  this  program. 
If  NDIM  or  any  NN(l)  is  less  than  one,  the  program  returns  immediately. 

Algorithm 

A  heavily  modified  version  of  the  algorithm  discovered  independently  by 
Danielson  and  Lanczos  ([2],  1942),  Good  ([4],  1958),  and  Cooley  and  Tukey  ([1], 
1965)  is  used.  The  following  example  is  an  application  to  a  one-dimensional 
transform  of  length  six. 

Let  w  =  e  The  transformation  is  written 

to  =  d-O  +  <*1  +  d.2  +  &3  +  <3.4  +  d.5 

ti  =  do  +  wdi  +  w2d.2  +  w3d3  +  w4cU  +  w5ds 
t>2.  =  do  +  v^di  +  w4d^  +  w6d3  +  w^d4  +  ‘W^'°d5 
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t3  =  do  +  w3di  +  w6da  +  w9d3  +  w12d  +  w15ds 

4 

t4  =  d<>  +  w4di  +  wsd2  +  w12d3  +  ■w16d4  +  w2°d5 
ts  =  do  +  v^dj.  +  wJ-°d2  +  w15d3  +  ■w2°d4  +  w25d5 

Straightforward  computation  requires  25  complex  multiplications  and  30  complex 
additions.  The  fast  Fourier  transform  computes  as  follows: 

uo  =  do  +  d3 
ui  =  do  +  w3d3 
u2  =  di  +  d4 
u3  =  di  +  w3d4 
u4  =  d2  +  ds 
U5  =  d2  +  w3d5 
to  =  Uo  +  U2  +  U4 
tl  =  Ul  +  WU3  +  W2Us 
t2  =  Uo  +  w 2u2  +  W4U4 

t3  =  Ux  +  \P\13  +  W6U5 

t4  =  Uo  +  w4u2  +  w8^ 

t5  =  Ul  +  \Pu3  +  wlou5 

which  requires  only  13  complex  multiplications  and  18  complex  additions.  Note 
that  w3  =  -1  and  w6  =  1. 

Such  a  reduction  in  computation  can  be  found  for  any  length  which  is  a 
composite  integer.  The  algebraic  proof  may  be  found  in  the  appendix.  Also, 
the  various  techniques  for  performing  multidimensional  transforms,  real  trans¬ 
forms,  etc.  are  discussed  there. 

Special  Cautions  and  Features 

The  finite  discrete  Fourier  transform  places  three  restrictions  upon 
the  data: 

1.  The  data  must  form  one  cycle  of  a  periodic  function.  Alternately 
stated,  the  subscripts  are  interpreted  modulo  N. 

2.  The  number  of  input  data  and  the  number  of  transform  values  must 
be  the  same. 
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3.  The  data  must  be  equispaced  in  each  dimension  (though,  of  course,  the 
interval  need  not  be  the  same  for  each  dimension).  Further,  if  in 
any  dimension  the  input  data  are  spaced  at  interval  dt,  the  resulting  transform 
values  will  be  spaced  from  0  to  2jt(N-1)/( Ndt)  at  interval  2n/(Ndt)  as  I  runs 
from  1  to  N.  By  periodicity,  the  upper  limit  is  identified  with  -2jt/(Ndt)  and 
in  fact  all  points  above  the  "foldover  frequency"  rt/(Ndt)  are  to  be  identified 
with  the  corresponding  negative  frequency. 

Those  familiar  with  other  implementations  of  the  fast  Fourier  transform 
may  be  aware  that  the  order  of  the  data  is  scrambled  in  the  course  of  execution. 
Unscrambling  is  performed  automatically,  however,  and  both  the  input  and  output 
values  are  placed  in  ordinary  sequential  arrangement. 

Timing 

Let  N,  ,  ,  be  the  total  number  of  points  in  the  data  array.  That  is, 
uOuaJ. 

Ntotal  =  N1*1*2***”  Decompose  into  its  prime  factors,  such  as 

K2  K'?  K5 

2  3  ?  ••••  Let  Zs  be  the  sum  of  all  the  factors  of  two  in  that  is, 

Zq  =  2*K2.  Let  be  the  sum  of  all  the  other  factors,  E_p  =  3*K3  +  5*K5  +  .... 
The  time  taken  for  a  multidimensional  transform  is 

T  =  To  +  Ntotal  [Tx  +  T2Ze  +  TfZf]  . 

For  the  CDC  3300, 

T  =  3000  +  N-total  [600  +  bOZe  +  175Lf]  microseconds. 

The  greater  optimization  apparent  for  factors  of  two  is  due  to 

1.  The  eight-fold  symmetry  of  the  trigonometric  functions  from  0  to  2«. 

2.  The  fact  that  Fourier  transforms  of  length  two  and  four  require 
fewer  complex  multiplies  than  transforms  of  other  lengths. 

The  above  timing  formula  is  accurate  for  complex  data. 

The  use  of  real  data  (IF0RM  =  0)  can  reduce  running  time  by  as  much  as 
forty  percent.  On  the  CDC  3300,  a  64  x  64  complex  array  was  transformed  in 
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6.1  seconds;  a  64  X  64  real  array  took  4.2  seconds.  A  complex  array  1500  long 
took  6.1  seconds,  while  a  real  1500  array  ran  only  3-4  seconds. 

Accuracy 

The  simplistic  idea  about  accuracy  is  apparently  correct:  because  the 
fast  Fourier  transform  takes  fewer  steps  in  execution,  less  error  creeps  in. 
Gentleman  and  Sande  ( [  3] ,  1966)  show  theoretically  that  the  root -mean- square 
relative  error  is  bounded  by 

1-06  4otal  2'b  V2fj]3/2 

where  b  is  the  number  of  bits  in  the  floating-point  fraction  and  f .  are  the 
*  0 
factors  of  Ntotal* 

Further  error  is  introduced  in  this  particular  program  by  the  use  of 
recursive  generation  of  sines  and  cosines  for  factors  of  Ntotal  other  than 
two.  Sines  and  cosines  needed  for  factors  of  two  are  computed  precisely.  In 
actual  practice,  out  of  eleven  and  a  half  digits  representable  on  the  CDC  3300, 
about  four  were  lost  on  long  one -dimensional  sequences  like  1500  and  4096. 

Applications 

Besides  all  the  direct  uses  of  discrete  Fourier  transforms  in  signal 
processing,  lens  design,  crystallography,  seismic  studies,  etc.,  Fourier 
transforms  find  application  in  techniques  of  correlation  and  convolution.  The 
principal  tool  here  is  the  convolution  theorem.  Denoting  the  convolution  of 
two  discrete  functions  f  and  g  by  f*g 

(f*g4  =  5  fA-3  * 

where  both  j  and  k  run  from  1  to  N  and  subscripts  are  interpreted  modulo  N, 
and  denoting  the  discrete  Fourier  transform  of  f  by  F(f),  the  convolution 
theorem  states 

F(f*g)  =  F(f)  F(g). 
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The  difficulties  here  are  that  cyclical  interpretation  of  subscripts  may 
not  be  desirable  and  that  N  may  not  be  convenient  for  fastest  processing  via 
the  fast  Fourier  transform.  Appendage  of  zeroes  to  the  ends  of  the  sequences 
solves  both  problems.  See  Stockham  ([5]>  1966)  and  Gentleman  and  Sande  ([3], 

1966) . 

Examples  of  Use 

A.  FOURT 

1.  Forward  transform  of  conqplex  50  X  40  array  in  Fortran  II 

DIMENSION  DATA  (2,50,40),  WORK  (100),  NN  (2) 

NN  (1)  =  50 
NN  (2)  =  40 
DO  1  1=1,  50 

DO  1  J  =  1,  40 
DATA  (l,I,j)  =  real  part 
1  DATA  (2,I,J)  =  imaginary  part 
CALL  FOURT  (DATA, NN, 2, -1,1, WORK) 

2.  Same  example  as  1,  but  in  Fortran  IV 

DIMENSION  DATA  (50,40),  WORK  (100),  NN  (2) 

COMPLEX  DATA 
DATA  NN/50,  40/ 

DO  1  I  =  1,  50 
DO  1  J  =  1,  40 
1  DATA  (l,j)  =  complex  value 

CALL  FOURT  (DATA,NN,2,-1,1,W0RK) 

3.  Same  example  as  2,  but  in  double  precision 

Add  the  following  statement: 

DOUBLE  PRECISION  DATA,  WORK 
Change  the  call  to: 

CALL  DFOURT  ( DATA, NN, 2, -1,1, WORK) 
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4.  Inverse  transform  of  real  64  x  32  array  in  Fortran  IV 

DIMENSION  DATA  (64,32),  NN(2) 

COMPLEX  DATA 
DATA  NN/64,32/ 

DO  1  I  =  1,  64 
DO  1  J  =  1,  32 
1  DATA(l,j)  =  real  value 

CALL  FOURT  (DATA, NN, 2, +1,0,0) 

B.  F0UR2 

Inverse  transform  of  real  64  x  32  array  in  Fortran  IV 

DIMENSION  DATA  (64,32),  NN(2) 

COMPLEX  DATA 
DATA  NN/64,32/ 

DO  1  I  =  1,  64 
DO  1  J  =  1,  32 
1  DATA(l,j)  =  real  value 
CALL  F0UR2  (DATA,NN,2,+l) 

C.  F0UR1 

Forward  transform  of  real  array  of  length  2048  in  Fortran  II 
DIMENSION  DATA  (2, 2048) 

DO  1  1=1,  2048 

DATA(l,l)  =  real  part 
1  DATA(2,I)  =  0 

CALL  F0UR1  ( DATA, 2048, -l) 
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Appendix  I 


Historical  Sketch 


In  1903  Runge  published  schemes  for  the  optimal  computation  of  twelve  and 
twenty-four  point  Fourier  transforms  ([6]).  They  involved  grouping  and  re¬ 
grouping  of  values  in  a  manner  similar  to  the  modern  FFT.  Runge 's  schemes 
are  well  known  and  appear  in  many  works  on  numerical  analysis,  including 
Runge  and  Konig  ( [7]  >  1924)  and  Whittaker  and  Robinson  ([8],  1944).  Neverthe¬ 
less,  no  one  thought  of  generalizing  Runge' s  ideas  until  1942  when  Danielson 
and  Lanczos  ([2])  published  an  optimal  algorithm  for  N  •  2^  point  transforms. 
Their  paper  passed  unnoticed. 

Meanwhile,  in  1937  Yates  ([9])  had  devised  an  algorithm  for  the  efficient 
computation  of  the  interactions  of  2n  factorial  experiments.  This  involves 

sums  of  the  form  .  .  .  .  . 

t.  =  Z  d.(-l)loJo+  llJl+“* 

J  1 

where  ioii  •••  and  joJi***  are  the  binary  representations  of  i  and  j. 

Davies  et  al  extended  the  method  to  3n  experiments  ([10],  1954);  three  years 
later,  Good,  in  an  abstruse  paper,  extended  it  to  general  factorial  experiments 
([4],  1958).  In  the  same  paper,  Good  devised  analogous  algorithms  for  N  point 
Fourier  transforms,  where  N  is  decomposable  into  mutually  prime  factors. 

Cooley  and  Tukey  removed  this  restriction  and  clarified  Good's  argument  ([1], 
1965).  They  also  wrote  what  was  probably  the  first  computer  program  to 
perform  FFT. 

Cooley  and  Tukey' s  paper  sparked  a  resurgence  of  interest  in  the  Fourier 
transform.  Despite  its  indispensability  in  many  areas  of  signal  processing, 
the  Fourier  transform  had  long  been  avoided  for  reasons  of  long  computation 
time.  The  FFT  revived  interest  to  such  an  extent  that  the  IEEE  Audio  Trans¬ 
actions  has  devoted  an  entire  issue  to  it  (June  1967)  and  three  groups  have 
proposed  implementing  it  in  hardware  ([11],  1963 ;  [12],  1967;  [13] ,  1967)* 
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Appendix  II 


The  Mathematics  of  the  Fast  Fourier  Transform 

Mathematical  descriptions  of  the  algorithms  used  in  the  Fast  Fourier 
Transform  subroutines  will  be  published  in  the  near  future. 

Punched  decks  for  these  three  subroutines  are  available  from  J.  J. 
Fitzgerald,  J-105,  or  from  SHARE. 
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Appendix  III 


Listing  of  the  Fortran  Subroutines 

The  listings  of  the  three  subroutines  F0UR1,  F0UR2,  and  FOURT  are  given 
on  the  following  pages.  All  three  are  written  in  USASI  Basic  Fortran,  and, 
as  such  are  compatible  with  the  great  majority  of  Fortran  compilers. 
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SUBROUTINE  FOuR1<DATA*NNfHION> 

""THE  OOOLBY'TUkSY  FAST  FOURIER  tRANEFORM 
TRANSFORM(J)  ■  SUM (DATA! I )*W**{ ( I*S)*( J»l) )  J  <  WHERE  !  AND  j  RUN 
FROM  1  TO  NN  AND  W  •  fiXR(lEJON*2*PI*SQRTUl)/NN).  DATA  IS  A  ONE- 
DIMENSIONAL  COMPLEX  ARRAY  <1*6,*  TWE  REAL  AND  IMAGINARY  RANTS  OF~ 
TH8  DATA  ARE  LOCATED  IMMEDIATELY  ADJACENT  JN  8T0RA0E*  IUCM  AS 
FORTRAN  IV  PLACES  THEM)  WHOSE  LENGTH  NN  IS  A  POWER  OF  TWO,  tSlQN 
IS  *1  OR  «t*  GIVING  THE  SION  OF  THE  TRANSFORM!  TRANSFORM  VALUES 
ARE  returned  IN  ARRAY  data*  replacing  THE  INPUT  DATA.  THE  TIME  It 
PROPORTIONAL  to  N*L0Q2(N)«  RATHER  THAN  THE  USUAL  N*«2,  WRITTEN  BY 
NORMAN  BRENNER,  jUNEtSST.  THIS  IS  THE  SHORTEST  VERSION 
OF  FFT  KNOWN  TO  THE  AUTHOR*  AND  IS  INTENDED  MAINLY  FOR 
DEMONSTRATION.  PROGRAMS  FOURS  AND  TOURT  ARE  -AVAILABLE  THAT  RUN 
-TRICE  AS  FAST  AND  OPERATE  ON  MULTIDIMENSIONAL  ARRAYS  WHOSE 
DIMENSIONS  ARE  NOT  RESTRICTED  TO  'POWERS  OF  TWO.  (LOOKING  UP  SINE9 
AND  COSINES  IN  A  TABLE  WILL  OUT  RUNNINO  TIME  OF  FOUR1  'EY  A  THIRD.) 
SIB--  IEEE  AUDIO  TRANSACTION®  SJUNi  SSSM* 

DIMENSION  DATaU) 

N»I*NN 

m 


!F(M-2)5* 3*3 
J«J*M 

JliAW 


THITA»J,t4199265SS*FL0AT(l8IQN*(M-l))/F.LOAT(MMAX) 

HR*COS(TH|TA> 


TIMPI«WR*DATA(J*l)*WI*DATAtJ) 
iDATA(J^»DATA<J>*T8MPR 


GO  TO  * 

'RETURN 

lNfr== 
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THE  COOUIV-TUKEY  r*JT  FOURIER  TRANEFOPM  IN  USUI  BASIC  FORTRAN 

TRANSPORT  Jl#  J2# , , , )  «  SUNIDATAI  U#  ft* # , ,J*W|**I  <  tiU>*  CJl-l) } 

WHERE  11  AND  Jl  RUN  PROM  1  TO  NNtll  AND  WHEXfrlt*f9N*S*PI* 
SQRT(«1)/NN(D if  ETCi 


C 

X 


fC‘ 

c 

■c 

m 

c 

X 

.'tr 

c 

c 

c 

c 

c 

c 

c 

c 


BATA~IS  A  HULT  IDlRlMf  16#Ai'  __  _  __ 

DIMENSIONS  ARE  ROWERS  OP  TWO,  tHi  .LENOTH  Of.jiACH  DIMENSION  I* 
STORED  IN  THE  INTEGER  ARRAY  NN,  OP  .I.8N0TH  NDIM,  IltON  T 
*1  OR  -1#  OIVINQ  THE  -SIGN  OP  THE 'TRANSFORM «  “THE  USAL 
AND  fMAQlNARY  PARTS  OP  A  iDATUH  ARE  IMMEDlATEt,?  ADJACENT  IN  8T0RA0E 
(SUCH  AS  FORTRAN  |V  PLACES  THEM),  TRANSFORM  IRSSUITS  ARE  RETURNED 
IN  ARRAY  DATA,  REPLACING  THE  ORIOINAL  DATA,  TIME  IS  PROPORTIONAL 
TO  N*LOQS(N)i  RATHER  THAN  THE  USUAL  NM3,  NOTE  THAT  IP  A  'FORWARD 
TRANSFORM  IS  FOLLOWED  EY  AN  INVERSE  TRANSFORM#  THE  ORDINAL  DATA 
KILL  REAPPEAR  MULTIPLIED  EY  NN U»*NNt2 >* , , , ,  EXAMPLE- 
FORWARD  FOURIER  TRANSFORM  OF  ,A  TWO-DIMENSIONAL  ARRAY  IN  FORTRAN  II 
DIMENSION  0ATA<2#64«32)iNN(2) 

NN(1)*64 


NN(2>»32 
DO  1  I *1#  64 
DO  1  J*l«  32 


DATA<1. Ii J>«REAL  PART 
0ATA<2# I|J>*IMAGINARY  PART 


C  CALL  P0UR21D*T*iKN»2»»l) 

C 

C  SAME  EXAMPLE  IN  FORTRAN  IV 

C  DIMENSION  DATA(64,32)#NNI2) 
c  complex  data 

C  DATA  NN/64# 32/ 

C  '  DO  1  l»l#64 

C  DO  1  J*li 32 

C  1  DATA ( I # J) iCOMpLEX  VALUE 
C  "CALL  F0UR2(DATA,NN#2#-1) 


C  PROGRAM  SY  NORMAN  BRENNER  FROM  THE  BASIC  PROGRAM  BY  CHARLES 
C  RADER,  may  1967,  THE  IDEA  FOP  THE  DIGIT  REVERSAL  KAS  SUGGESTED 
C  BY  RALPH  ALTER, 

C  ‘THIS  VERSION  OP  ThE  f AST  FOURIER  TRANSFORM  IS  THE  FASTEST  -KNOWN 
C  TO  THE  AUTHOR,  LOOKING  UP  SINES  , AND  COSINES  IN  A  TABLE  INSTEAD  OF 
C  COMPUTING  THEM  WOULD  DECREASE  RUNNING  TIME  SEVEN  PERCENT, 

C  PROGRAMS  FOURT  AND  FOUR1  ARE  AVAILABLE  FROM  THE  AUTHOR  THAT  ALSO 

C  PERFORM  THE  PAST  FOURIER  TRANSFORM  AND  ARE  WRITTEN  IN  USAS!  BASIC 

C  FORTRAN,  FOURT  IS  THREE  TJMES  AS  LONG#  IS  NOT  RESTRICTED  TO 

C  POWERS  OF  TWO,  AND  RUNS  UP  TO  FORTY  PERCENT  FASTER  ON  REAL  DATA, 

C  FOUR1  IS  ONE  FOURTH  AS  LONG#  ONE  HALF  AS  FAST#  AND  JS  RESTRICTED 

C  TO  ONE  DIMENSION  AND  POWERS  OF  TWO,  _ 

C  SEE--  IEEE  AUDIO  TRANSACTIONS  (JUNE  19671#  SPECIAL  ISSUE  ON  FFT, 

C 

DIMENSION  DATA(1)«NN(1> 

IF(NO!M*>1)700,1#1 
1  NT0T«2 

DO  2  IDIH«1#NDIM 
IF(NNUDIli)  )700#700#2 
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r2  NTOT*NTOf*NN<lDJM> 

RTHLF*, 70710  67612 
Twopufyassie  53070 

C  MAIN  LOOP  /OR  each  dimension 

c 

NP1»2 

DO  600  IDIM«1»NDIM 
N«NN< I D I M ) 

NP2«NPl*N 

IF ( N-l >  700 . 600 . 100 
C 

C  SHUFFLE  DATA  0Y  BIT  REVERSAL!  SINCE  N»2**K,  AS  THE  SHUFFLING 
C  CAN  BE  DONE  By  SIMPLE  INTERCHANGEi  NO  WORKING  ARRAY  IS  NEEDED 
C 

100  NP2HF»NP2/2 


J»1 

DO  160  1 2»1, Np2| NP1 
IFU-I2)110*130,130- 
110.  IiMAX«I2*NPl-2 

DO  120  li»I2i I1MAX.2 
DO  120  1 3«  11# nTOT, NP2 
J3* J* 1 3”  1 2 
TEMPR»DATA<  J3) 
TEMPl«DATA< 13*1) 

D AT A  < 1 3 ) ■DAT  A{ J3 ) 

DATA(I3*l)*DATA(J3*l> 
DATA(J3)*TEMPR 
120  DATA  I J3*l ) *TEmP I 
130  M»NP2HF 
140  IFU*M>160»160*150 

150  g«j-M 
M»M/2 

IFIM-NPU160.140.140 


160  J«J*M 
C 

c  main  loop,  perform  fqurjer  transforms  of  length  four,  with  one  of 

C  LENGTH  two  IF  NEEDED,  THE  TWIDDLE  FACTOR  W«GXPl !S1GN*2*PI* 

0  SQRT(-1)*M/(4*MMAX)),  CHECK  FOR  THE  SPECIAL  CASE  W»JS!QN«SQRT{*1> 

C  AND  REPEAT  FOR  W«W*H*lSlGN*SQRTt»l>  )/SQRT(2) , 

C 


NP1TW»NP1*NP1 

TPAR’N 

310  IFIIPAR«2)350, 33 0,320 
320  IPAR»IPAR/4 
GO  TO  310 


330 


DO  3<0  llil.Npl.2 
DO  340  K1«H,NT0T,NP1TW 
*K2»K1*NP1 
•tsmPR»DATA(K2> 
TEMP1»DATA(K2*1) 
DATA<K2)«DATA(Kt).TEMPR 


DATA<K2*U«DATA(Kl*l)«TEMPJ 
DATA<Kl)»DATA(Kl)*TEMPR 
340  ^DATAI Kl*l )*D *T A(K1*1 J*TEMP1 

350  MMAX»NPl 

360  IFIMMAX»NP2HF)3?0, 600,600 
LMAX»M*X0 (NPtTW.MMAX/2) 

DO  570  LrNP1«LMAX,NP1TW 
M»L 


W!«S!N(TMITA) 

•410  W2MWft*HR-Wl*Hl 

us  i  *a , 

w3MM2R«WRpW8I*W! 

W3!»W8R*WI*N2t*K* 

Tljfl - DO  5 JO  !lal«Npi'2 

KMtN»lRAMMMl 
f(MM4¥8NPl>430#430»440 
N«U  -  ■  ~~  " 

f«tPAR«HHAX 

Rl»KHJN*NTOT#KlTiP 
KI«Kl*KDir 
Kl*Kl*KOir 


*4t0 


U«MDATA<«1  _ 
U81»DATA(KJ*1>*DATAJK4*1) 


U4|«DATA(K4»*DATA(KJ) 

00  TO  *10 

OAlfl 
Oitft 

8ll«H2R*0ATA<K8)a8t1«DATA(Ra«l) 
Ttti«8A*BATA(Ka*l>*N8l*0ATA(K2> 
TJMWR*0ArA(K3)«Wi«DATA(NJ*l) 
*KM0*fUKJU»*kf!*6iTA<l<l> 

_ _ Ti*IE»i 

fi(7A(Kl)*T2ft  == 
Uli*0ATA(KlAl>*T8| 

U8R«TJA*T4R 


'U4MTJHT4I 

U4t«T4R«T3R 

BQ  iv 


,fOO  y4iif 4  t  gf 3 1 

WWW  Win  T  T I  f  T  ■  I 


rAtKipiMutt^uai 

rA(Kl)*uaR»U4R 


ISO 


»  'AtK4*M*U?t«iU4| 
N»4M4HINMA>*U 
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SUBROUTINE  rOURT1piTSTNNTNPfHT|ttSNiffORHrWORK) - — 

THE  C001.IY-TUKEY  FAST  FOURIER  TRANIFORM  IN  USAS!  BASIC  FORTRAN 

TRANSFORM* J1#J2# , , , )  «  8UM(-DATA<  U#  12#  , , , >*W1**( ( ll«l >*<Ul-i) ! 

WHERE  U  AND  Jl  &UN  FROM  1  TO  NNti)  ANO  Wl«EXfr(tStGN*t*Pt» 
SQRTU1)/NNU)>«  ETC,  THERE  IS  NO  .UNIT  ’ON  THE  DIMENSIONALITY 
(NUMBER  OF  SUBSCRIPTS!  OF  THI  DATA  ARRAY,  IF  AN  INVERSE 
TRANSFORM  USjQNMl)  IS  PERFORMED  (IRON  AN  ARRAY  OF  TRANIFORMBD 
(lS!ONa*l)  DATA,  THE  ORIOINAL  DATA  HILL  REAPPEAR. 

MULTIPLIED  BY  NN(1)*NN(2)*.|,  THE  ARRAY  OF  INPUT  DATA  MUST  IE 
IN  COMPLEX  FORMAT,  HOWEVER.  IF  Alt  tMAOlNARY  PARTS  ARB  2BR0  (!,f, 
THE  DATA  ARE  DlSOUJSBD  REAL)  RUNNING  'T  tME  IS  CUT  UP  TO  FORTY  PER¬ 
CENT,  (for  fastest  transform  of  .real  data,  nnui  should  be  even.) 

THE  TRANSFORM  VALUES  ARE  ALWAYS  COMPLEX,  ANO  ARE  RETURNED  IN  TH* 
ORIGINAL  ARRAY  ;OF  DATA,  REPLACING  THE  INPUT  DATA,  THE  LENSfH 
OF  EACH  DIMENSION  OF  THE  DATA  ARRAY  MAY  BE  ANY  INTEGER,  THE 
PROGRAM  RUNS  FASTER  ON  COMPOSITE  INTBOIRS  THAN  ON  PRIMES,  AND  IS 
PARTICULARLY  FAST  ON  NUMBERS  .RICH  IN  FACTORS  OF  TWO, 

"TIMING  IS  IN  FACT  GIVEN  IT  THE  FOLLOWING  FORNUL*#  Lt?  NTOT  IE  VHE 
TOTAL  NUMBER  OF  POINTS  (REAL  OR  COMPLEX)  IN  THE  DATA  .ARRAY,  THAt 
IS,  NT0T*NNU)*NN<2)*,.,  DiOOMPOSE  NTOt  INTO  ITS  PRIMS  FACTORS, 
SUCH  AS  2**K2  *  3**K3  *  S**K5  ♦  ,,,  .LIT  SUN2  BE  THE  SUM  OF  ALL 
THE  FACTORS  OF  TWO  IN  NTOT,  THAT  IB,  SUMS  ■  2*K2,  LIT  9UMF  Si 
THE  SUM  OF  ALL  OThBR  FACTORS  >0F  NTOT#  THAT  IS#  SUHF  ■  3*W3*»*K5* , , 
THE  TIME  TAKEN  BY  A  MULTIDIMENSIONAL 'TRANSFORM  ON  THESE  ’NTOT  'DATA 
IS  T  «  TO  ♦  NT0T*(Ti*T2*8UM2*T3*SUMF)|  ON  THE  OOC  3300  (FLOATING 
POINT  AOD  TIME  ■  SIX  MICROSECONDS)#  T  •  JOflO  ’*  NTOT*( *00*4Q«SUm3* 
173*SUMF>  MICROSECONDS  ON  COMPLEX  DATA, 

IMPLEMENTATION  OF  THE  DEFINITION  :8Y  SUMMATION  HILL  RUN  IN  A  TIME 
PROPORTIONAL  TO  NT0T*(NNU)*NNt2 >♦,,.) ,  FOR  HIGHLY  ‘COMROSlTl  NTOT 
THE  SAVINGS  OFFERED  BY  THIS  PROGRAM  CAN  BE  DRAMATIC*  A 'ONB«DImEN> 
SIONAL  ARRAY  4000  {N  LENGTH  WILL  BE  TRANSFORMED  IN  4000*(60Q* 
40*(2*2*2*2*2)*173*(5*5*5)>  «  14,5  SBCONDS  i'ERSUS  ABOUT  4000* 
4000*175  *  2800  SECONDS  FOR  THE’ STRAIGHTFORWARD  TECHNIQUE , 

THE  FAST  FOURIER  TRANSFORM  ■’PLAGES  THREE  RESTRICTIONS  UPON  THE 
DATA. 

1,  the  number  of  input  data  and  the  number  of  transform  values 
must  be  the  same, 

2,  BOTH  the  input  data  and  the  transform  values  must  rbpresent 
eouispaced  points  in  their  respective  domains  of  time  .and 
frequency,  calling  these  spacings  oeltat  and  deltaf#  it  must  be 
true  THAT  DELTAF*2*Pl/(NN( I )*DILTAT) ,  OF  COURSE.  1DELTAT  NEED  NOT 
BE  THE  SAME  FOR  EVBRY  DIMENSION. 

3,  conceptually  at  least#  the  Input  data  and  the  transform  output 

REPRESENT  SINGLE  cycles  OF  PERIODIC  FUNCTIONS. 

THE  CALLING  SEQUENCE  is»- 

CALL  FOURT(DATA,NN#NDIM. {SIGN# IFORM.WORK) 

DATA  is  THe  ARRAY  USED  TO  HOLD  THE  REAL  AND  IMAGINARY  -PARTS 
OF  the  data  On  INPUT  AND  THE  TRANSFORM  VALUES  ON  OUTPUT,  IT 
IS  A  MULTIDIMENSIONAL  FLOATING  POINT  array#  WITH  THE  REAL  AND 
IMAGINARY  PARTS  OF  A  DATUM  STORED  IMMEDIATELY  ADJACENT  (N  STORAGE 
(SUCH  AS  FORTRAN  IV  PLACES  THEM),  ’NORMAL  FORTRAN  ’ORDER ING  IS 
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ncinoooonoonooano 


c  expected*  the  FIRST  subscript  changing  pastes?,  the  dimensions 

C  ARE  GIVEN  IN  THE  INTEGER  ARRAY  NN,  OP  LENGTH  NDIM.  ISJQN  JS--1 
C  TO  INDICATE  A  FORWARD  TRANSFORM  (EXPONENTIAL  SIGN  IS  •)  AND  •*! 

C  FOR  AN  INVERSE  TRANSFORM  (SIGN  IS  *),  IPORM  IS  *1  IP^THE  DATA  ARE 

C  COMPLEX*  0  IF  THE  DATA  ARE  REAL,  IF  IT  IS  0,  THE  JMAO JNARY 

C  PARTS  OF  THE  DATA  MUST  BE  SET  TO  .ZERO)  AS  EXPLAINED  ABOVE,  THE 

C  TRANSFORM  VALUES  ARE  ALWAYS  COMPLEX  AND  ARE  STORED  IN  ARRAY  DATA, 

C  WORK  IS  An  ARRAY  USED  FOR  WORKING  STORAGE,  IT  IS  FLOATING  POINT 

C  REAL*  ONE  DIMENSIONAL  OF  LENGTH  EQUAL  TO  TWICE  THE  LARGGS?  ARRAY 

*C  DIMENSION  NN ( I )  THAT  IS  NOT  A  POWER  OF 'TWO,  IF  ALL  4NN< I )  ARE 

C  POWERS  OF  TWO,  IT  JS  NOT  NEEDED  AND  MAY  BE  REPLACED  BY  ?lRO  IN  THE 

c  calling  sequence,  thus,  for  a  one-dimensional  array,  nn<i>  odd, 

C  WORK  OCCUPIES  AS  MANY  STORAGE  LOCATIONS  AS  DATA,  IF  SUPPLIED* 

C  WORK  MUST  NOT  BE  THE  SAME  ARRAY  AS  -DATA,  ALL  SUBSCRIPTS  OP  ALL 

C  ARRAYS  BEGIN  AT  ONS, 

-e— 

c  example  i.  three-dimensional  forward  fourier  transform  >of  a 

C  COMPLEX  ARRAY  DIMENSIONED  32  BY  25  BY  13  I N  FORTR'AN  IV. 

C  DIMENSION  DATA(32,25»13)*HORK(EO)iNN(3) 

c  v  complex  data 
c  DATA  NN/32,29,13/ 

C  DO  1  1*1*32 

C  DO  1  J*1 • 25 

C  DO  1  KU,13 

C  1  DATA<1.JiK)*COMPLEX  VALUE 
C  -CALL  F0URT(DATA,NN,3r-l,l,W0RK) 

C  ' 

C  EXAMPLE  2,  OnE-DJMENS JONAL  FORWARD  TRANSFORM  OF  A  REAL"tARRAY  OF 

C  LENGTH  64  IN  FORTRAN  II, 

C  DIMENSION  DAT  A( 2, 64 ) 

C  DO  2  I *1*  64 
C  DATA(1* I )«R8*L  P*RT 
C  2  DATA(2*l>«0, 

C  CALL  FOURT<D*T*,64,1,-1,8,0> 

C 

C  THERE  ARE  NO  ERROR  MESSAGES  OR  ERROR  iHALTS  IN  THIS  PROGRAM,  ThE 
program  returns  immediately  IF  NDIM  OR  ANY  NN C I i  IS  LESS  than  ONE, 

PROGRAM  BY  NORMAN  BRENNER  FROM  THE  BASIC  PROGRAM  BY  CHARLES 

RADER,  JUNE  1967.  THE  IDEA  FOR  THE  'DIGIT  REVERSAL  WAS 
SUGGESTED  BY  RALPH  ALTER, 


THIS  IS  THE  FASTEST  AND  MOST  VERSATILE  VERSION  OF  THE  FFT  KNOWN 
TO  THE  AUTHOR,  A  PROGRAM  CALLED  F0UR2  IS  AVAILABLE  THAT  A^SO 
PERFORMS  THE  FAST  FOURIER  TRANSFORM  AND  IS  WRITTEN  IN  USA$I  BASIC 
FORTRAN,  IT  IS  ABOUT  ONE  THIRD  AS  LONO  AND  RESTRICTS  THE 
DIMENSIONS  OF  T^E  INPUT  ARRAY  (WHJOH  MUST  BE  COMPLEX)  TO  BE  POWERS 
OF  TWO,  ANOTHER  PROGRAM,  CALLED  FQUR1,  IS  ONE  TENTH  AS  LONG  AND 
RUNS  TWO  THIRDS  AS  FAST  ON  A  ONE-DjMBNSlONAL  COMPLEX  ARRAY  WHOSE 
LENGTH  IS  A  POWER  OF  TWO, 


REFERENCE-- 

IEEE  AUDIO  TRANSACTIONS  (JUNE  1967),  SPECIAL  ISSUE  ON  THE  FFT, 

DIMENSION  DAT A(l),NNli)f IFACT( 32) *WORKlD 
TW0Pl»6.a831853Q7 


RTHLF*, 78716  67812 
t F<ND IMn 1)92 0,1*  i— 
NT0T»2 

DO  2  IDIH«1,NDIM 
IF (NN( I  DIM) )920 ,920, 2 
NTOT»NTOT*NN(lDIM>  . 
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— 


TC“ 

c 


HinW  UO0P  FOR  EACH  DIMENSION 


NP1*2 

DO  «0  III 
N*NN< i D I M  > 
NP2»NRltN 

"lffifiimiivBv** 


c 

c 


IS  N  A  POWER  or  TWO  AND  IF  NOT,  WHAT  ARE  ITS  FACTORS 


iPil 

!DIV»2 

10  IOUOT«H/lD!V 

f OCMVM  ' 

In*"  n 


ifactufuidiv 

IF* IF*1 


[NON* • 

!OUOT*H/|DIV 
IRIH*H*ID I V*lOUOT 


SI 

11 


miipotp  tB  i  v  mrnm 

IF(IREMM0»38f4tr? 

!PACTUPl*IDIf§ES 

ir»rr*t  ^ 

M* IOUOT 
QO  TO  30 


91 


irt|RiH>*0»91»A0 

NTW0»NTW0*NTW0 

<00  TO  70 

IFACTtlFAtH 


MOOltfliHjOfltil** 
tT 

DIMENSIONS, 

2,  RIAL  TRANSFORM  FOR  THE  2ND  DR  3RD  DlMBNllON, 


_ 


IT 


IMirnlN 
340  tMi'rrM/trAomn 

— ■ 
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•850  jij-irw 

rrp2«tr»l  .  ^•^-•--^-  ■----  -■  ■ 

if*if*i 

imFp8»NP17260i?6ffi*«0 


860 


870 

C 

C 

C 

C 

C 

C 

300 

305 

310 

320 

330 


CONTINUE 

!2MAXal3*NP2-NPl 
1*1 

DO  270  tl«  J3|  I2MAXiNP1 
DATA  I J2)  *H0RK(  I ) 

DATA(r8+»)«H0RK(l*l) 

Ht*2 

MAIN  COOP  FOR  FACTORS  OF  TWO,  PERFORM  FOURIER  TRANSFORMS  OF 
LENQTH  FOUR*  WITH  ONE  OF.CINSTM  TWO  IF  NEEDED*  THE  ‘TWIDDLE  FACTOR 
W«BXP(ISIQN*8*PI*S0RT(.U*M/(4*MMA¥>J1  OHECOOR  W» JSION*$QRTt-l> 
AND  REPEAT  FOR  W«W*{t*i$I0N*80RT(-lJt/S0RTt2)1 

IF(NTWO«NP1)600, 600*309 
NPJTW»NP1*NP1 

ipar»ntwo/npi 

IF<IPAR-8)350, 

IPAR«JPAR/4 
OO  TO  310 

DO  340  U«1*URNQ«8 
DO  340  Kl»ll»NTOT*NPlTW 


340 

350 

360 

370 


K8»K1*NP1 

TEMPR»0ATA<K2) 

’TEMPI»DATaIK2*1) 

DATA(K2J*DATA(K1)-TEMPR 

DATA<K2*1>»DATA(K1*1>«TEMPI 

DATAUi.)«DATA(K17*TEMPR 

DATA(K1*1>»DATA(K1<i1>*TEHP! 

MMAX*NP1 

IF (MMAX*NTWO/2 1370*600,600 
LMAX«MAX0(NPiTW,MMAX/2) 

DO  570  L*NP1#LMAX,NP1TW 


M*U 

IF(MMAX«NP1J420*420,380 
380  THKT Ai-TWOP  i*FLOAT (L> /FLOAT! 4*MMAX) 
IFCISIGNMOO, 390*390 
390  THETAi-THETA 
400  WRaCOS(THETA) 

WMSlNtTHETA) 

410  W2R»WR*WR»Hl*Wl 
W2I,2,*WR*WI 
W3R»W2R*WR-W2l*WI 
W3I»W2R*Wl*W8l*WR 
420  DO  530  llal, llRNQ*8 
KM1N«I1*IPAR*M 
IF(MHAXpNP1)430»430,440 
430  KMIN*U 
440  *KDIFsIPAR*MM*X 
450  KSTEP»4*»KDIF 

IF!KSTEPpNTWO)460,460*530 
460  DO  520  Kl«KMlN#NTOTiKSTEP 
K2aKl*KDlF 
K3sK2*KDIF 
K4»K3*KDIF 

IF<MMAX^NP1)470, 470,480 
470  U1R»0ATA«K1)*DATA(K2) 

U1!»DATA(K1*17+DATA<K2*1> 

U2R»DATA<K37*DATA(K4) 
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<71 

<72 

<80 


<90 

800 

510 


520 


530 


540 

550 


560 


570 


C 

•e 

C 

c 

c 

X 

600 


U2I«DaTA<k3*1 I<DATA<K4<U 
U3R«DATA<K1>’DATA<K2> 

U3I*DATA<K1*1)-DATA(K2*1) 

IFC ISI0N><71i472»472 
U4R»DATA<K3*1)-0ATA<K4*1) 

U4I«DATA<K<)-D*TA(K3) 

00  TO  510 

U<R»DATA<K<*U"DATA(K3*1) 

U< J«DATA(K3)-DATA(K<) 

00  TO  510 

T2R»W2R<DATA<k2>-W2 J<DATA<K2<U 
T2!«W2R*DATA<k2*1)*W2I*DATA<K2) 

T3R<WR*DAT A <K3>«W J<DATA<K3<U 
T3I<WR<DATA<K3<U<WI<0ATA  <KJ> 

T4R*W3R*DATA<K<)”W3I <DATA<K4<U 
T4t <W3R<DATA<K<<U <H3l<DATA?K<) 

U1R<DATA<KU<T2R 
U1I<DATA<K1<U<T21 
U2R«TJR<T«R 

U21<T31<T<I 
U3R<DATA<KU<T2R 
U3t«DATA<Kl*U-T2l 
!F( ISIGN)490*500»500 
U4R»T3t«T4l 
*U4I"T4RsT3R 
00  TO  310 
U<R<T<HT3I 
U4t*T3R*T4R 
DATA<KU<U1R<U2R 
DATA<K1<U<U1  J<U2l 
*DATA<K2><U3R<t|4R 
DATA<K2<U<U3J<U<1 
DATA<K3)«U1R*U2R 
DATA<K3<U<U1J<U2I 
DATA<K<)»U3R-U6R 
DATA<K4<U<U3J<U4J 
KDJF»KST6P 
KMlN«4*<KHIN»lU*ll 
00  TO  <50 
CONTINUE 

m»m*lmax 

!F(M-MMAX>340, 540,370 
trUS10N>550«360»360 
TEMPR<WR 

WR<<WR<WII<RTHLF 
WI<<WJ-TEMPR>*RTHLF 
00  TO  .410 
TEMPR<WR 

*r.<WR-WU*RTmL* 

wi«<tehpr*wi)«rthlf 

00  TO  410 
CONTINUE 
JPAR<3<IPAR 
MMAX<MMAJ<<MMAX 
00  T®  360 

MAIN  LOOP  FQR  FACTORS  NOT 

W<BXP<ISIQN<2<PI<S0RT<<U  _  ...  _ .....  .  ... 

PERFORM  A  FOURIER  TRANSFORM  OF  .LBNOTH  IFACTUFJi  MARINO  lUSB  OF 
CONJUGATE  SYMMETRIES, 

tf(NTWO*Mp2}605|7Q0i7Q0 
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6ft  IF«*lfTI8- 


610 


615 


620 

625 


lf»TVDI« 

NPlHF«NPl/2 
rrP2MPACT(  jF)*irpt 
-j£fUNtaPi«| 


690 


635 

640 

645 

650 


ir<JlMJN»!FPi>6l5,6l5«640 
00  639  Jl*diMlNj IFPlrNPl 
TRETA»»TW0Pt,Fl0mJi»l>/ru0AT<tPP2) 
tr(ISIQN)625<620i6t0 
THBTA»'THET4 
»5TPR«COS{THETAJ 
WSTP!«SIN<TH6T*> 

WRsWSTPR 
«l«WSTPI 
j2MiN»ji*irpi 
J2MAX«ji*!FP2.irPl 


PS  ___ 

IlMAX«J2*llPNa*9 
DO  630  il#j2<  JlMA 
DO  630  j3«Ji«NT0T, 

TBHP"»DATA(J3> 

DATA  t  j3)«D*TA(  J3)*WR«DATA(  J3*i1>*N1 


TEMPP-WH 

WR«MR*USTPR»tt|*HSTPI 
Vl«TEMPR*H5TPi*WI*WSTPR 
THBTA.-TWOPl/rtOATUFACTUr) ) 
tf(JSlGN) 69 0*645*645 
TRfTA.-THETA 


USTPRaCOS(THSTA) 

WSTP!»S JN(ThETAI 
J2RNO«irPi*U*irACTUP)/I> 
DO  695  11*1* llRNO, 2 
DO  695  (S«U«NTOT|NP2 
J2HAX4I3*J2RN0»  JTP1 
00  690  J3«|3*J2MAX*IPP1 
JlMAX-J2*IfPl-NPl 

DO  600  Jl»j2* jlMAX»NPl 
JJMAX«J1*NP2*!FiP2 
DO  600  J3bJ1* J3HAX* IPP2 
JHJN*J39J2*I3 
JMAX«JMlN*|FP2»|PPl 
I«**{J3e|3>/NPiHr 
IF (J2*t 3)659 *695 *665 
9UHR*0 *  * 

SUMl'Ot 

DO  660  JPAX* |FP1 

SUMR"8UMR*DATaCj) 
SUMl450MJ*DAtitJ*t) 


655 


>460 


1>I0RK(1)«SUMR 


665 


WORK  t t*l)»8UMJ 
00  TO  610 


JfJRAX 

suhr*datau» 

OLDSR«0  J 
OLDS J«0 1 

«*f*r 


670 


PHI®  .. 

Tf HPRiSUMR 
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169  WRaWSTPP 

_ B 

W|«TINPR*MSTPl*Wl*WSTPR 
690  TWOHR«WR*WR 
iiij 

§  ttBSXBlfGSRlAMRl. 

HO  695  {*•!}» j2H*X»NPi 

DATA(t2*l)*NORK(I«l) 

695  !■!*• 

irpi«*rP2 

C  COMPLETE  A  real  TRANSFORM  in  THE  .1ST  'DIMENSION#  N  SVEN#  'BY  CON- 
C  'JUGATE  SYMMETRIES, 

700  GO  TO  (9O0»8OO#9OO»7OX)» ICASS 

701  NMALFpN  =-  -  — 

- Nf  - - 

tmsta»»twopi^fioat(n>- 

irUStQN>703«7Qa#702  • 

‘702  THiTAp-TMETA  - 

•703  wSTPR«CO»(TMETA) 

NSTPIf Sf N{TH6T6I 
Nt»RSTPt 
WlANSTPJ 
IMIN*3 

I  fJMTN»t*N#ALF’ 

■*80  TO  72f 
710  fJljMIN 

•M"7*0  IflMfNiNTOTif 
9UMR»<DATa<  l  UDATaUH/2, 

SUMl"(DATA(I*l>*DATA(jPl))/2, 


iDATA«I)?SUMR*tEMPR 
DATA«r*l)fDirt»TlMRl 


D*T*tJ*llA'-OIFI*TEMPI 
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'7*8  g»J*N»a 


-1MIN»1MIN»2 

JM!N«JH!N"2 

TEMPRbWR 

WRbWR«wSTpRbWI*H8TPI 
W|«TEHPR*WSTPI*¥I*¥STPR 
725  tFUMlN«jM!N>7lQi730,74Q 
738  IF<ISIQN>731»74Q,740 

73i  00  735  T*!MiN#NT0T4WP 

735  DATA<tn>*»DATA(**U 
740  NP2»NP2*NP2 

ntot»ntot*ntot 

J«NTOT*l 


745 


758 

755 


|MAX*NT0T/2*1 
!M!N»IMAX-2*NHALF 
I«IMlN 
QO  TO  755 
BATA<J)!BATA<J? 

OATAl J+l)«*DATA(|*U 
HI+8 
g«j-2 


ir < I-IMAKJ7S0# 760,740 
760  DATA(J>9DATA(:RIN>-DATAUM!N*n 
DATAtJ*l?»0, 

m  ira-j)77o*7Bo#78o 

765  DATAlJ)«OATAU> 

DATA* J*^7  «DAT *t 1*1 > 

770  j • t«a 

g»j*a 


DATA< J*l>«8, 
IMAX*IMjN 
00  TO  745 


780  OATA ||>|8 ATA (1 UPATA (2 ) 

COHPUETG  A  real  TRANSFORM  FOR  THE  2ND  OR  3RD  DIMENSION  BY 
e  CONJUGATE  SYMMETRIES, 

0 

^ftj=IFmRNO«NPt>609,900»»08 
80S  DO  860  U«1iNT0T*NP2 
!2HAX«I3*NR2'>NR1 

—  DO  860  raBl3,l2MAX,NRi== 

IMtN<l2*URN0 

IMAX»l2*Npl«2 

JHAX"2*r3*NPl-tM!N  - 

IFna»I3)820i  820,818^^ 

110  jMAXBjHiX#NP2 

820  TfUOtMil)850, 850,830  —— 
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